/**
 * @file test_smft.cpp
 *
 * Test of SFMT
 *
 * @author Dahua Lin
 */

#include <gmlpp/core/random.h>
#include <gmlpp/core/accum.h>

#include <gmlpp/base/timer.h>
#include <cstdio>

using namespace gml;

void calc_stat(const colvec& x)
{
	std::printf("\t mean = %.4f, var = %.4f\n", mean(x), var(x));
}

void dd_stat(uint32_t K, const Col<uint32_t>& x)
{
	Row<uword> counts(K);
	counts.zeros();

	intcount(srange(K), x, counts);

	std::printf("p = ");
	for (uint32_t k = 0; k < K; ++k)
	{
		std::printf("%.3f ", double(counts[k]) / x.n_elem);
	}
	std::printf("\n");
}


void test1(default_rng_t& rg)
{
	const uword nrep = 10;
	const uword N = 1000000;

	colvec x(N);

	// uniform distribution

	std::printf("Uniform distribution:\n");

	Timer tm0;
	tm0.start();

	for (uword k = 0; k < nrep; ++k)
	{
		fill_randu(rg, x);
	}

	double et0 = tm0.elapsed_usecs();
	std::printf("\t throughput = %.4f MPS\n", N * nrep / et0);
	calc_stat(x);

	// normal distribution

	std::printf("Normal distribution:\n");

	Timer tm1;
	tm1.start();

	for (uword k = 0; k < nrep; ++k)
	{
		fill_randn(rg, x);
	}

	double et1 = tm1.elapsed_usecs();
	std::printf("\t throughput = %.4f MPS\n", N * nrep / et1);
	calc_stat(x);

	// exponential distribution

	std::printf("Exponential distribution:\n");

	Timer tm2;
	tm2.start();

	for (uword k = 0; k < nrep; ++k)
	{
		fill_rande(rg, x);
	}

	double et2 = tm2.elapsed_usecs();
	std::printf("\t throughput = %.4f MPS\n", N * nrep / et2);
	calc_stat(x);

	// gamma distribution

	std::printf("Gamma distribution (alpha = 2):\n");

	Timer tm3;
	tm3.start();

	for (uword k = 0; k < nrep; ++k)
	{
		fill_randg(rg, x, 2.0);
	}

	double et3 = tm3.elapsed_usecs();
	std::printf("\t throughput = %.4f MPS\n", N * nrep / et3);
	calc_stat(x);

}


void test2(default_rng_t& rg)
{
	const uword nrep = 10;
	const uword N = 1000000;
	const uint32_t K = 4;
	double p[K] = {0.3, 0.2, 0.4, 0.1};
	double w[K] = {3, 2, 4, 1};

	Col<uint32_t> x(N);

	std::printf("ddsample:\n");

	Timer tm0;
	tm0.start();

	for (uword k = 0; k < nrep; ++k)
	{
		for (uword i = 0; i < N; ++i)
			x[i] = ddsample(rg, K, p);
	}

	double et = tm0.elapsed_usecs();

	std::printf("generation speed = %.4f MPS\n", N * nrep / et);
	dd_stat(K, x);

	std::printf("ddsample_w:\n");

	Timer tm1;
	tm1.start();

	for (uword k = 0; k < nrep; ++k)
	{
		for (uword i = 0; i < N; ++i)
			x[i] = ddsample_w(rg, K, w);
	}

	double et1 = tm1.elapsed_usecs();

	std::printf("generation speed = %.4f MPS\n", N * nrep / et1);
	dd_stat(K, x);
}


void test3(default_rng_t& rg)
{
	const uword nrep = 10;
	const uword N = 1000000;
	WORSampler<uword> sp(N);

	Col<uword> x(N);

	Timer tm1;
	tm1.start();

	for (uword k = 0; k < nrep; ++k)
	{
		for (uword i = 0; i < N; ++i)
		{
			x[i] = sp.next(rg);
		}
		sp.reset();
	}

	double et = tm1.elapsed_usecs();

	std::printf("wor sample speed = %.4f MPS\n", N * nrep / et);

	Col<uword> y(N);
	for (uword i = 0; i < N; ++i) y[i] = i;

	Timer tm2;
	tm2.start();

	for (uword k = 0; k < nrep; ++k)
	{
		rand_suffle(rg, y);
	}

	double et2 = tm2.elapsed_usecs();

	std::printf("suffle speed = %.4f MPS\n", N * nrep / et2);
}


int main(int argc, char *argv[])
{
	default_rstream_t rstream;
	default_rng_t rg(rstream);

	test1(rg);

	return 0;
}



